Global definitions

Probability-generating functions

t0=Sys.time()
#Definition of functions

h <- function(z) {
  (1 - z)^((1 - z) / z)
}

px <- function(z,Am,m1,Amm) {
  (h(z)^{Am}) * (h(h(z)^m1)^{Amm})
}

psb <- function(z,pb) {
  (1 - pb) + pb * z
}

ps <- function(z,p) {
  (1 - p) + p * z
}

pgb <- function(z,Am,m1,pb,Amm) {
  px(psb(z,pb),Am,m1,Amm)
}

pgs <- function(z,p,Am,m1,Amm){
  px(ps(z,p),Am,m1,Amm)
}
#Requiered library for 3D plots
library(plotly)

#Print twenty decimal places
options("digits" = 20)

#Data set to work with
FILENAME="R-GM - FinalPop 2^33, bottle neck 2^20, growth cycles 3, WT2R 10^-7, WT2M 1e-06 M2R 5e-05.txt"

#FILENAME="R-QM - FinalPop 2^33, bottle neck 2^20, growth cycles 3, WT2R 10^-7, WT2M 1e-05 M2R 5e-04.txt"
#FILENAME <- "luria_4trans_1.00e-007mua_2seed.txt"
#Reading the file
td <- read.table(FILENAME,header = T,sep="")

#Sampling probability
ps1 <- 2^{20}/2^{33}

#Dilution process: Binomial sampling, with probability ps1, of resistant cell in the j-th growth cycle
binomial_sampling_GC_j<-function(j,ps1){
  a=c()
  for (i in 1:nrow(td)){
    a=c(a,rbinom(n = 1, size = (td[, (4 * (j - 1) + 2)]+td[ ,(4 * (j - 1) + 4)])[i], prob = ps1))
  }
  return(a)
}

#Recovery of the k-th probability state of the pmf
#WT2R  Mutation rate from wildtype to resistant cells
#WT2M  Mutation rate from wildtype to mutator cells
#pop   Final population size
#N     Total number of probability states
pmf_coeff<-function(k,WT2R,WT2M,pop,N){
  sum=0
  for(n in 0:(N-1)){
    #Auxiliary variable to transform the pgf into a characteristic function
    aux_char=exp((2i*pi*n)/N)
    #Auxiliary variable to apply the FFT algorithm
    aux=(-2i*pi*n*k)/N
    #FFT                                   #Fixed scaling factor: 500
    sum=sum + pgs(aux_char,p=ps1,Am=10^{WT2R}*pop,m1=10^{WT2M},Amm=10^{WT2R}*pop*500)*exp(aux)
  }
  #Real part
  return(Re(sum/N))
}

First growth cycle

#Binomial sampling of resistant cells in the j-th growth cycle, j=1
data<-binomial_sampling_GC_j(j=1,ps1)

#Maximum value in data 
N=max(data)
N
## [1] 33
steps=30

#log-likelihood function
ll_fun<-function(WT2R,WT2M){
  sum(log(pmf_coeff(data, WT2R,WT2M,pop=2^33,N)))
}

paramters_WT2R=c(-9,-5)
paramters_WT2M=c(-8,-4)

output=c()
for(i in 0:steps){
  for(j in 0:steps){
    m_WT2R=paramters_WT2R[1]*(paramters_WT2R[2]/paramters_WT2R[1])^{i/steps}
    m_WT2M=paramters_WT2M[1]*(paramters_WT2M[2]/paramters_WT2M[1])^{j/steps}
    ll=ll_fun(m_WT2R,m_WT2M)
    output=rbind(output, c(10^{m_WT2R}, 10^{m_WT2M}, 500*10^{m_WT2R}, ll))
  }
}
colnames(output)=c('WT2R', 'WT2M', 'M2R', 'log-likelihood')
head(output)
##                           WT2R                      WT2M
## [1,] 1.0000000000000000622e-09 1.0000000000000000209e-08
## [2,] 1.0000000000000000622e-09 1.5230713629662131529e-08
## [3,] 1.0000000000000000622e-09 2.2975616244345304987e-08
## [4,] 1.0000000000000000622e-09 3.4334922203825658472e-08
## [5,] 1.0000000000000000622e-09 5.0841703260677619694e-08
## [6,] 1.0000000000000000622e-09 7.4612270297717860124e-08
##                            M2R         log-likelihood
## [1,] 5.0000000000000008327e-07 -744.49903306033786521
## [2,] 5.0000000000000008327e-07 -744.49115671835045305
## [3,] 5.0000000000000008327e-07 -744.47966875932934272
## [4,] 5.0000000000000008327e-07 -744.46306921362315734
## [5,] 5.0000000000000008327e-07 -744.43930306396680407
## [6,] 5.0000000000000008327e-07 -744.40558319719457359
#Maximum
maximum_GC1=output[order(output[,4], decreasing=TRUE),][1,]
maximum_GC1
##                       WT2R                       WT2M 
##  7.6856069542669374823e-08  8.0631442982151470284e-05 
##                        M2R             log-likelihood 
##  3.8428034771334688524e-05 -2.6653791929043785559e+02

Log-likelihood surface

#options(width = 10000, height=1000 )
discretization<-function(range, n){
  mu_min=range[1]
  mu_max=range[2]
  a=seq(0,n)
  result=mu_min*(mu_max/mu_min)^{a/n}
  return(result)
}

paramters_WT2R=c(-9,-5)
WT2R=discretization(paramters_WT2R,steps)
paramters_WT2M=c(-8,-4)
WT2M=discretization(paramters_WT2M,steps)

zz=matrix(output[,4],nrow = steps+1,ncol=steps+1)

x.lab="Mutation rate from WT to R on log<sub>10</sub> scale"
y.lab="Mutation rate from WT to M on log<sub>10</sub> scale"
z.lab="log-likelihood"
main=""


x=-7
y=-6
z=ll_fun(-7,-6)
df_real<-data.frame(x, y, z)

x=log10(maximum_GC1[1])
y=log10(maximum_GC1[2])
z=maximum_GC1[4]
df <- data.frame(x, y, z)


#plot_ly(showscale = T, name=" ", showlegend = TRUE)

fig=plot_ly(x=~WT2R, y=~WT2M, z=~zz, type = 'surface',name= "Log-likelihood surface", showlegend = TRUE, colorbar = list(title = "Log-likelihood"), contours = list(z = list(show=F, usecolormap=TRUE,highlightcolor="#ff0000",project=list(z=T)))) %>%  layout(autosize = F, width = 900, height = 650,
      title = paste("\n",main),
      scene = list(
        xaxis = list(title = x.lab),
        yaxis = list(title = y.lab),
        zaxis = list(title = z.lab)
      ))%>%add_trace(data = df, x = x, y = y, z = z, mode = "markers", type = "scatter3d",name= "Maximum", showlegend=T ,marker = list(size = 6, color = "red", symbol = 104))

fig=fig%>%add_trace(data = df_real, x = x, y = y, z = z, mode = "markers", type = "scatter3d",name= "Real values", showlegend=T ,marker = list(size = 5, color = "deepskyblue", symbol = 104))

fig
df=as.data.frame(output)


df_confidence<-df[df$`log-likelihood` >= maximum_GC1[4]-qchisq(0.95,2)/2, ]

min(df_confidence$WT2M)
## [1] 1.0000000000000000209e-08
max(df_confidence$WT2M)
## [1] 0.00010000000000000000479
fig1=fig%>%add_surface(showscale = F,name="95% confidence region" ,
  contours = list(
    z = list(
      show = TRUE,
       project=list(z=T),       # (don't) project contour lines to underlying plane
       usecolormap = F,        # (don't) use surface color scale for contours
      color = "red",             # set contour color
      width = 30,                   # set contour thickness
      highlightcolor = "black",  # highlight contour on hover
      start = min(df_confidence$`log-likelihood`),                   # include contours from z = 0...
      end = max(df_confidence$`log-likelihood`),                  # to z = 1400...
      size = max(df_confidence$`log-likelihood`)-min(df$`log-likelihood`)                   # every 100 units

    )
  )
)
fig1=fig1%>%layout(autosize = F, width = 900, height = 650)
fig1
#fig %>% add_trace(z = matrix(min(df$`log-likelihood`), nrow = nrow(zz), ncol = ncol(zz)), contours = list(coloring = "none",  color = "black"), showscale = F, name="Lower limit of the 95% confidence region")
 
#fig

#fig%>%add_trace(x=~log10(df$WT2R), y=~log10(df$WT2M), z=~df$`log-likelihood`, type='scatter3d',  mode = "markers", marker = list(size = 2))

#fig


# Calculate the posterior probabilities
#df$posterior_prob <- (df$`log-likelihood`)  / sum((df$`log-likelihood`))

# Sort the data frame by posterior probability in descending order
#data_frame <- df[order(-df$posterior_prob), ]

# Calculate the cumulative posterior probability
#data_frame$cumulative_prob <- cumsum(data_frame$posterior_prob)

# Identify the values within the 95% credible region
#credible_region <- data_frame[data_frame$cumulative_prob <= 0.975, ]
#credible_region <- credible_region[credible_region$cumulative_prob >= 0.025,]




#fig%>%add_trace(x=~log10(credible_region$WT2R), y=~log10(credible_region$WT2M), z=~credible_region$`log-likelihood`, type='scatter3d',  mode = "markers", marker = list(size = 2))

#fig

Likelihood surface

#options(width = 10000, height=1000 )
discretization<-function(range, n){
  mu_min=range[1]
  mu_max=range[2]
  a=seq(0,n)
  result=mu_min*(mu_max/mu_min)^{a/n}
  return(result)
}

paramters_WT2R=c(-9,-5)
WT2R=discretization(paramters_WT2R,steps)
paramters_WT2M=c(-8,-4)
WT2M=discretization(paramters_WT2M,steps)

zz=matrix(exp(output[,4]-max(output[,4])),nrow = steps+1,ncol=steps+1)

x.lab="Mutation rate from WT to R on log<sub>10</sub> scale"
y.lab="Mutation rate from WT to M on log<sub>10</sub> scale"
z.lab="Rescaled likelihood"
main=""

x=log10(maximum_GC1[1])
y=log10(maximum_GC1[2])
z=exp(maximum_GC1[4]-max(output[,4]))
df <- data.frame(x, y, z)


#fig=plot_ly(showscale = T, name=" ", showlegend = TRUE)

fig=plot_ly(x=~WT2R, y=~WT2M, z=~zz, type = 'surface',name= "Rescaled likelihood surface", showlegend = TRUE, colorbar = list(title = "Rescaled likelihood"), contours = list(z = list(show=F, usecolormap=TRUE,highlightcolor="#ff0000",project=list(z=T)))) %>%  layout(autosize = F, width = 900, height = 650,
      title = paste("\n",main),
      scene = list(
        xaxis = list(title = x.lab),
        yaxis = list(title = y.lab),
        zaxis = list(title = z.lab)
      ))%>%add_trace(data = df, x = x, y = y, z = z, mode = "markers", type = "scatter3d",name= "Maximum",
                     marker = list(size = 6, color = "red", symbol = 104))

fig=fig%>%add_trace(data = df_real, x = x, y = y, z = exp(z-z), mode = "markers", type = "scatter3d",name= "Real values", showlegend=T ,marker = list(size = 5, color = "deepskyblue", symbol = 104))


fig=fig%>%add_surface(showscale = F,name="95% confidence region" ,
  contours = list(
    z = list(
      show = TRUE,
       project=list(z=T),       # (don't) project contour lines to underlying plane
       usecolormap = F,        # (don't) use surface color scale for contours
      color = "red",             # set contour color
      width = 20,                   # set contour thickness
      highlightcolor = "black",  # highlight contour on hover
      start = exp(min(df_confidence$`log-likelihood`)-max(df_confidence$`log-likelihood`)),                   # include contours from z = 0...
      end = exp(max(df_confidence$`log-likelihood`)- max(df_confidence$`log-likelihood`)),                  # to z = 1400...
      size = exp(max(df_confidence$`log-likelihood`)- max(df_confidence$`log-likelihood`)) -exp(min(df_confidence$`log-likelihood`)-max(df_confidence$`log-likelihood`))                   # every size units
    )
  )
)
fig=fig%>%layout(autosize = F, width = 1100, height = 650)
fig